Hugoniot of shocked liquid deuterium up to 300 GPa: Quantum molecular dynamic simulations 
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Quantum molecular dynamic (QMD) simulations are introduced to study the thermophysical properties 
of liquid deuterium under shock compression. The principal Hugoniot is determined from the equa- 
tion of states, where contributions from molecular dissociation and atomic ionization are also added 
onto the QMD data. At pressures below 100 GPa, our results show that the local maximum compres- 
sion ratio of 4.5 can be achieved at 40 GPa, which is in good agreement with magnetically driven flyer 
and convergent-explosive experiments; At the pressure between 100 and 300 GPa, the compression ra- 
tio reaches a maximum of 4.95, which agrees well with recent high power laser-driven experiments. In 
addition, the nonmetal-metal transition and optical properties are also discussed. 
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Recent studies of materials under extreme conditions, 
which require improved understandings of the thermophys- 
ical properties in the new and complex regions, have gained 
much scientific interest yj]. The combination of high tem- 
perature and high density defines "warm dense matter" 
(WDM) - a strongly correlated state, which is characterized 
by partially dissociated, ionized and degenerated states, 
and the modelling of the dynamical, electronic, and opti- 
cal properties for such system is rather challenging. Due 
to their simplicity, hydrogen and its isotopes (deuterium 
and tritium) have been studied intensively [2], where the 
relative pressure and temperature have reached megabar 
range and several eV. Specially, as one of the target ma- 
terials in the inertial confinement fusion experiments J3], 
deuterium has been extensively investigated through exper- 
imental measurements and theoretical models. Gas gun 
[4], converging explosive Oa], magnetically driven flyer 
[6], and high power laser-driven experiments have 

been applied to probe the physical properties of deuterium 
during single or multiple dynamic compression. Theoret- 
ically, approximations have been introduced to simulate 
warm dense deuterium (hydrogen), such as, linear mixing 
model [10], chemical model FVT [11], path integral Monte 
CarloJPIMC) and quantum molecular dynamics 

[Hlill- 
To date, although a number of explanatory and predictive 
results in some cases have already been provided by experi- 
mental and theoretical studies, however, many fundamental 
questions of deuterium under extreme conditions are still 
yet to be clarified. The equation of states (EOS), especially 
the Hugoniot curve, are essential in this context. Since five 
to six-fold the initial densities have been detected bylaser- 
driven experiments at megabar pres sure regime OHM] and 
supported by PIMC simulations Ill2lll4ll . considerable con- 
troversies in the deuterium EOS have been raised. Mean- 
while, converging explosives [5] and magnetically driven 
flyer ]6J] experiments indicate that the compression ratio 
(r/) shows a maximum close to 4.3, which is in good agree- 
ment with QMD results 1 15, H]. Furthermore, in adiabatic 



and isentropic compressions, the nonmetal-metal transi- 
tion of deuterium (hydrogen), which is accompanied by the 
change of optical spectroscopies [17], has been a major is- 
sue recently. The links between nonmetal-metal transition 
and dissociation (ionization) under dynamic compression 
are of particular significance yfl. 

The chemical pictures of deuterium under extreme con- 
ditions could be briefly described as two processes: (i) par- 
tial dissociation of molecules, D 2 <^ 2D, and (ii) a subse- 
quent ionization of atoms, D^ e+D + . QMD simulations, 
where electrons are modelled by quantum theory, are con- 
vinced to be a powerful tool to describe the chemical reac- 
tions, such as dissociation and recombination of molecules. 
Meanwhile, the dynamical, electrical and optical proper- 
ties of warm dense matter have already been proved to 
be successfully investigated by QMD simulations IU8U19I1 . 
However, the ionization of atoms is not well defined in the 
framework of density functional theory (DFT). Consider- 
ing these facts, thus, in this paper we applied the corrected 
QMD simulations to shock compressed deuterium, and the 
calculated compression ratio is substantially increased ac- 
cording to the ionization of atoms in the warm dense fluid. 

We have performed simulations for deuterium by em- 
ploying the Vienna Ab-initio Simulation Package (VASP) 
|2(A ulll . A fixed volume supercell of N atoms, which is 
repeated periodically throughout the space, forms the ele- 
ments of the calculation. By involving Born-Oppenheimer 
approximation, electrons are fully quantum mechanically 
treated through plane -wave, finite-temperature DFT II 1511 . 
and the electronic states are populated according to the 
Fermi-Dirac distribution at temperature T e . The exchange 
correlation functional is determined by generalized gra- 
dient approximation (GGA) with the parametrization of 
Perdew-Wang 91 I22I1 . The ion-electron interactions are 
represented by a projector augmented wave (PAW) pseu- 
dopotential [23]. The system is calculated with the isoki- 
netic ensemble (NVT), where the ionic temperature T t is 
kept constant every time step by velocity scaling, and the 
system is kept in local thermodynamical equilibrium by 
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setting the electron (T e ) and ion (Tj) temperatures to be 
equal. 




FIG. 1 : (Color online) Calculated pair correlation function (black 
line) and coordination number (red line) at temperatures of 3000 
K (solid line), 5000 K (dashed line), and 10000 K (dotted line). 
Inset is the contour plot of the ionization degree as a function of 
density and temperature. 

The plane -wave cutoff energy is selected to be 600.0 eV, 
so that the pressure is converged within 5% accuracy. T 
point is employed to sample the Brillouin zone in molecu- 
lar dynamic simulations, because EOS (conductivity) can 
only be modified within 5% (15%) for the selection of 
higher number of k points. A total number of 128 atoms 
(64 deuterium molecules) is included in a cubic cell, and 
over 300 (densities and temperatures) points are calculated. 
The densities adopted in our simulations range from 0.167 
to 0.9 g/cm 3 and temperatures between 20 and 50000 K, 
which highlight the regime of principal Hugoniot. All the 
dynamic simulations are lasted for 4 ~ 6 ps, and the time 
steps for the integrations of atomic motion are 0.5 ~ 2 fs 
according to different densities (temperatures). Then, the 
subsequent 1 ps simulations are used to calculate EOS as 
running averages. 



TABLE I: Coefficients in expansion for the internal energy 
E. 
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In QMD simulations, zero point vibration energy 
(■^his V ib) and ionization energy (13.6 eV/atom) are ex- 
cluded, thus, the internal energy and pressure should be 
corrected as follows: 

1 



E = E, 



qmd + g JV(1 - a)E mh + NfiEi, 



P = Pqmd + (1 + P) 



pk B T 



(1) 



(2) 



where N is the total number of atoms for the present sys- 
tem, and m D presents the mass of deuterium atom. The 
density and temperature are denoted by p and T respec- 
tively, and k B stands for Boltzmann constant, a and j3 are 
the dissociation degree and ionization degree, while E vib 
and E ion correspond to the zero point vibration energy 
and ionization energy, respectively. Eq M d and Pqm d are 
calculated from VASP. Various corrections to QMD sim- 
ulations have already been applied to model warm dense 
matter lfl5l [l6Tl . but contributions from atomic ionization, 
which are particularly important at high pressure, are still 
in absence. Ionization degree of aluminium under ex- 
treme conditions has been successfully quantified through 
Drude model II24I1 . however, the simple metallic model 
is not suitable for the present system. Here, a new and 
effective method in accounting for contributions to the 
EOS from molecular dissociation and atomic ionization 
has been demonstrated. 



TABLE II: Coefficients b ik in expansion for the total pressure P. 
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The dissociation degree is important in determining the 
internal energy, from which EOS can be derived, especially 
at low temperatures and the initial state on the Hugoniot 
curve. The dissociation degree could be evaluated through 
the coordination number: 
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where O, is the volume of the supercell. The coordina- 
tion number is a weighted integral over the pair correlation 
function (PCF) g(r) of the ions. The doubled value of K at 
the maximum of g(r) (r = 0.75 A), is equal to the fraction 
of atoms forming molecules in the supercell. The sampled 
PCF and K(r) are labelled in Fig. Q] Fast dissociation of 
molecules emerges at the temperature between 4000 and 
7000 K, and a region featured with (dP/dT) v < 0, which 
is not presented here, is observed. Our results show that 
molecular deuterium can be neglected above 15000 K due 
to thermal dissociation. 

The ionization degree of the system can be evaluated 
through Saha equation: 
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where only one level of ionization process is considered. 
In the present formula, m e stands for the electron mass. 
As shown in the inset in Fig. [TJ the ionization of deu- 
terium could be neglected below 10000 K. At the temper- 
ature between 10000 and 50000 K, where the modeling of 
the principal Hugoniot of deuterium is rather difficult, par- 
tially ionized warm dense fluid is formed, and the ioniza- 
tion of atoms is of predominance in determining the EOS. 
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FIG. 2: (Color online) Simulated principal Hugoniot curve (solid 
red curve). Previous data are also shown for comparison. The- 
ories: QMD results by Lenosky et al. fl^l (dotted blue curve) 
and Hoist et al. 1 16] (dashed blue curve); restricted PIMC sim- 
ulations by Magro et al. fl2ll (dashed black curve), Militzer et 
al. 1 1 3fl (dotted black curve), and direct PIMC by Bezkrovniy 
et al. Y\M (dashed dot black curve); the linear mixing model of 
Ross iflOTI (dashed green curve), and the chemical model FVT 
I Un (dotted green curve). Experiments: gas gun by Nellis et al. 
|4fl (solid circle), Z-pinch by Knudson et al. (g] (solid square), 
explosives of Boriskov et al. [5[] (open circle), laser-driven by 
Hicks et al. | 9] (up open triangle) and Boehly et al. (3] (down 
open triangle). 



Following Lenosky et al. 02511 . Beule et al. [26], and 
Hoist et al. lfl6Tl . we fit the internal energy and pressure 
by expansions in terms of density (g/cm 3 ) and tempera- 
ture (10 3 K). The corrected QMD data for internal energy 
(eV/atom) can be expanded as follows: 
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The total pressure given in GPa can be similarly expanded 
as E with the expansion coefficients bj k . The expansion 
coefficients a ik and bj k for E and P (accuracy better than Ci(cj) 
5%) are summarized in Tab. Hand Tab. HU respectively. 

Based on the EOS, the principal Hugoniot curve can be 
derived from the following equation: 



(E -E 1 ) = k---)(P + P 1 ), 
2 Po Pi 



where the subscripts and 1 refer to the initial and shocked 
states. In our present simulations, the initial density p 
is 0.167 g/cm 3 with the respective internal energy E = 
—3.28 eV/atom at T = 20 K. The pressure P of the 
starting point on the Hugoniot can be neglected compared 
to high pressures of shocked states. 

The principal Hugoniot is shown in Fig. [2] where previ- 
ous theoretical and experimental results are also provided 
for comparison. At pressures below 100 GPa, our results 
indicate that the principal Hugoniot experiences a local 
maximum compression ratio of 4.5 around 40 GPa, which 
can be attributed to the dissociation of molecules. The 
present Hugoniot agrees well with previous experiments, 
such as gas gun [4], magnetically launched flyer plates (6J], 
and converging explosives [5]. At the pressure between 40 
and 100 GPa, the Hugoniot curve shows a stiff behavior, 
and the compression ratio lies between 4.25 and 4.5. Mean- 
while, the ionization of atoms increases remarkably at P > 
50 GPa (see the inset in Fig. and consequently intener- 
ates the fluid. Thus, the combined effect of the molecular 
dissociation and atomic ionization results in a local mini- 
mum of r] (4.25) along the Hugoniot. 

Recent high power laser-driven experiments sug- 
gest that deuterium is stiff (rj max ~ 4.2) below 100 GPa 
and become softer (rj ~ 4.5 ~ 5.5) above 110 GPa, which 
can be described by the present simulations. Our results 
indicate that molecular deuterium can be neglected at this 
stage, and the atomic ionization dominates the character- 
istic of the Hugoniot with t] lies between 4.5 and 4.95 
(maximum is reached at 200 GPa), which is accordant 
with recent experiment [27]. The wide -range behavior of 
the Hugoniot is characterized by two stage transitions — 
dissociation under low pressure and ionization at higher 
pressure, and the present results show excellent agreement 
with experimental ones. The Hugoniots from mere QMD 
simulations hardly exceed 100 GPa I115ll except for that of 
Hoist et al. jlvH . but r] does not exceed 4.5 (at P > 100 
GPa). Although some PIMC simulations Oj, [H]] show 
five to six-fold compressions, the simulated data are not yet 
comparable with experiments. Due to the intrinsic approx- 
imations, no consistency has been detected between our re- 
sults and those of linear mixing model lllfJll and chemical 
model FVT ifTTTl. 

Let us turn now to see the nonmetal-metal transition by 
studying the optical and conductive behaviors of the warm 
dense deuterium. The real part of dynamic conductiv- 
ity CTi(w) can be evaluated through the following Kubo- 
Greenwood formula: 
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where i and j summations range over N discrete bands 
included in the calculation. The a sum is over the three 
spatial directions. /(e i; k) describes the occupation of the 
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FIG. 3: 



(Color online) Calculated dc conductivity along the 
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Hugoniot curve (solid squares). Previous data 0, ll^l are also 
shown for comparison. Inset is the ionization degree along the 
Hugoniot. 
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FIG. 4: (Color online) Calculated optical reflectivit y o f 
length 808 nm along the Hugoniot. Previous data fj ^fl7l . l28ll 
also plotted for comparison. 



and previous experiments lfl7ll . The increase of reflectance 
(from 0.05 to 0.6) is observed, and this can be interpreted 
as a gradual transition from a molecular insulating fluid to 
a partially dissociated and metallic fluid at above 40 GPa. 

In summary, we have performed QMD simulations to 
study the thermophysical properties of deuterium under ex- 
treme conditions. The Hugoniot EOS has been evaluated 
through QMD calculations and corrected by taking into ac- 
count the molecular dissociation described by the coordi- 
nation number K(r) and the atomic ionization described 
by Saha equation. The corrected Hugoniot has shown good 
agreement with the experimental data in a wide range of 
shock conditions, which thus indicates the importance of 
physical picture of a two-stage transition, i.e., dissocia- 
tion and ionization. The principal Hugoniot reveals a lo- 
cal maximum compression ratio of 4.5 at 40 GPa. With 
the increase of pressure, r] reaches a maximum of 4.95 
at about 200 GPa, and the contribution from atomic ion- 
ization demonstrates softened character of the Hugoniot. 
Smooth transition from a molecular insulating fluid to an 
partially dissociated and metallic fluid are observed at 40 
GPa. Our calculated optical constants along the Hugoniot 
have shown excellent agreement with experiments. In ad- 
dition, smooth fit functions constructed in the present paper 
for the internal energy and total pressure are expected to be 
useful for the future studies of warm dense deuterium. 

This work was supported by NIFC. 



ith band, with the corresponding energy e.^k an d the wave- 
function k at k. w(k) is the k-point weighting factor. 

From the calculated dynamic conductivity along the 
Hugoniot curve, the dc conductivity which follows 
the static limit oj— >0 of cri(u), is then extracted and plot- 
ted in Fig. [3] as a function of the Hugoniot pressure. As 
shown in Fig. [3] a^ c increases rapidly with pressure up 
to 40 GPa towards the formation of metallic state of deu- 
terium, which agrees well with the experimental measure- 
ments [4]. Similar tendency has also been found in the 
QMD simulations of warm dense hydrogen lUfll . When 
further increasing the Hugoniot pressure, one finds from 
Fig. [3]that a dc keeps almost invariant and the warm dense 
deuterium maintains its metallic behavior. Here we address 
that the nonmetal-metal transition is induced by gradual 
dissociation of molecules and thermal activation of elec- 
tronic states, instead of atomic ionization, which is not ob- 
served until 50 GPa according to the charge density distri- 
bution in the QMD simulations. Quantitative analysis can 
be clarified through plotting the ionization along the Hugo- 
niot as shown in the inset in Fig. [3] Meanwhile, optical 
reflectivity, with the respective wavelength of 808 nm, is 
shown along the principal Hugoniot in Fig. [4] where good 
agreement has been achieved between our present work 



* Corresponding author; zhang_ping@iapcm.ac.cn 

[1] R. Ernstorfer et al, Science, 323, 1033 (2009). 

[2] W. J. Nellis, Rep. Prog. Phys., 69 1479 (2006). 

[3] F. Philippe et al, Phys. Rev. Lett., 104 035004 (2010). 

[4] W. J. Nellis et al., J. Chem. Phys., 79 1480 (1983). 

[5] G. V. Boriskov et al, Phys. Rev. B, 71 092104 (2005). 

[6] M. D. Knudson et al., Phys. Rev. B, 69 144209 (2004). 

[7] G. W. Collins et al, Science, 281 1 178 (1998). 

[8] T. R. Boehly etal, Phys. Plasmas, 11 L49 (2004). 

[9] D. G. Hicks et al, Phys. Rev. B, 79 0141 12 (2009). 
[10] M. Ross, Phys. Rev. B, 58 669 (1998). 
[11] H. Juranek etal, J. Chem. Phys., 112 3780 (2000). 
[12] W. R. Magro et al, Phys. Rev. Lett., 76 1240 (1996). 
[13] B. Militzer et al, Phys. Rev. Lett., 85 1890 (2000). 
[14] V. Bezkrovniy et al, Phys. Rev. E, 70 057401 (2004). 
[15] T. Lenosky et al, Phys. Rev. B, 61 0163-1829 (2000). 
[16] B. Hoist et al, Phys. Rev. B, 77 184201 (2008). 
[17] P. M. Celliers etal, Phys. Rev. Lett., 84 5564 (2000). 
[18] A. Kietzmann et al, Phys. Rev. Lett., 101 070401 (2008). 
[19] W. Lorenzen etal, Phys. Rev. Lett., 102 115701 (2009). 
[20] G. Kresse et al, Phys. Rev. B, 47 R558 (1993). 
[21] G. Kresse et al, Phys. Rev. B, 54 1 1 169 (1996). 
[22] J. P. Perdew, Electronic Structure of Solids (Akademie Ver- 

lag, Berlin, 1991). 
[23] P. E. Blochl, Phys. Rev. B, 50 17953 (1994). 
[24] S. Mazevet etal, Phys. Rev. E, 71 016409 (2005). 
[25] T. J. Lenosky et al, Phys. Rev. B, 56 5164 (1997). 



5 



[26] D. Beule etal, Phys. Rev. B, 59 14177 (1999). [28] L. A. Collins etal, Phys. Rev. B, 63 184110 (2001). 

[27] M. D. Knudson et al, Phys. Rev. Lett., 103 225501 (2009). 



